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Model fields from the Los Alamos National Laboratory 
(LANL) Parallel Ocean Program (POP) 1/6° global circulation 
model are compared to measurements from over 1300 satellite- 
tracked surface drifters that were deployed in the tropical 
Pacific (20°N to 20°S) , between 1979 and 1994, during the 
TOGA Pan-Pacific Current Study. Geographic averages of 5-day 
averaged drifter velocity estimates for 2° latitude x 8° 
longitude bins are compared to similarly binned 3 -day model 
snapshots from September 1992 to October 1994 . 

Eulerian comparisons of the model mean velocities and 
their observed counterparts show that the model u mean is 
slightly higher in the equatorial region, while the model v 
mean is 50% greater in the this region. Model SST mean 
values are 20% less than observed values in the eastern 
equatorial Pacific. Model variability is about 20% less than 
the observed quantity in equatorial regions, and 50% less 
poleward of 10°S and 10°N. Both model and observed velocity 
and SST covariance fields imply a net heat convergence toward 
the equator- with the largest values in the region of 
instability waves north of the equator. 

Model velocity fields are used to produce simulated 
Lagrangian trajectories for uniform and non-uniform 
deployment strategies . Autocorrelation, time and length 
scales, diffusivity, and polarization are calculated and 
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ensemble averaged by 5° latitude bands for comparison with 
drifter-based Lagrangian statistics. Time and length scales 
are too long and dif fusivities are too low compared to 
observations, but data sampling in the simulated fields was 
biased by trajectories that overlap current regimes. 

These differences, in both Eulerian and Lagrangian 
comparisons, may be related to the lack of a surface mixed 
layer, inadequate representation of wind forcing, still too 
coarse grid resolution, and deficiencies in simulating the 
mean structure of the density field in the model. They are 
also partly related to lack of weighted averages to account 
for non-uniform drifter sampling. 
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INTRODUCTION 



I . 



High-resolution eddy-resolving models have shown 
increasing realism over the past decade, contributing to our 
understanding of ocean circulation and its variability 
(Hurlburt et al . , 1992; Semtner and Chervin, 1992; Smith et 
al . , 1992; Bryan et al . , 1995; Bleck et al., 1995; Semtner 
1997; Maltrud et al., 1997). Validation of these models is 
obtained by comparing the model output with high-quality data 
sets, such as those collected during the World Ocean 
Circulation Experiment (WOCE) and from satellites. Prior to, 
and during WOCE, reliable surface drifter technology in the 
Pacific Ocean has produced current and temperature 
measurements over spatial and temporal scales sufficient for 
statistical comparisons with high resolution models, 
providing a measure of the model's ability to reproduce the 
observed mean flow and variability (Niiler et al., 1997) . 

In the Pacific, surface drifter measurements have been 
collected by several investigators over the past two decades . 
In the early eighties, McNally (1981, 1983) deployed drifters 
in the mid-latitudes of the North Pacific, while Hansen and 
Paul (1984) used them to study currents associated with 
barotropic instability waves in the eastern equatorial 
Pacific . 

In 1987, as part of the Tropical Ocean and Global 
Atmosphere (TOGA) Program, surface drifters were deployed 
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across the entire Tropical Ocean basin to directly measure 
the circulation of the tropical Pacific, and to investigate 
temporal variability over scales from intraseasonal to 
interannual. This Pan-Pacific Current Study (Niiler et al . , 
1997) became part of the WOCE program objective to deploy 
several thousand drifters worldwide, providing a global 
surface data set of sea surface temperature, velocity, and 
atmospheric pressure. The study involved the deployment of 
over 1300 satellite-tracked surface drifting buoys, drogued 
to 15 m, in the tropical Pacific between January 1988 and 
December 1994, between 20°S and 20°N. The data set spans the 
1991-92-93 El Nino, the 1988 La Nina, and three normal 
seasonal cycles . This data set was used for both Eulerian 
and Lagrangian analyses, and has been characterized as 
appropriate for comparison with OGCM's of the tropical 
Pacific which have been forced with 1988-1994 atmospheric 
fluxes . 

Key studies of surface drifter tracks conducted from 
this data set, using Lagrangian statistical analyses, include 
a description of the seasonal variability of the equatorial 
Pacific, in Reverdin et al. (1994), and the interannual 
variability of the equatorial Pacific, in Frankignoul et al . 
(1996) . Studies of this data set using Eulerian analysis 
include Niiler et al., (1997), and estimates of horizontal 
divergence and vertical velocity in the equatorial Pacific, 
presented in Poulain (1993) . 
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The objective of this study is to calculate and compare 
the observed Eulerian and Lagrangian drifter statistics from 
the tropical Pacific with those from the Los Alamos National 
Laboratory (LANL) Parallel Ocean Program (POP) model. 

Eulerian averages in 2° latitude by 8° longitude sub-regions 
are compared directly, while a fourth-order Runge-Kutta 
scheme is used to generate trajectories from model velocity 
fields, from which Lagrangian statistics are calculated and 
related to their drifter-based counterparts. Comparison 
results are used to characterize inconsistencies for further 
model refinement. 

This study is organized as follows. An overview of the 
tropical Pacific is presented in Section II. Descriptions of 
the drifter data set and model output, along with background 
summaries of model and drifter investigations are presented 
in Section III. Drifter and model comparisons of time mean 
circulation, temperature fields, and variability are 
presented in Section IV. Section V describes temporal and 
spatial scales, measures of angular momentum and 
polarization, and dif fusivities as a result of Lagrangian 
statistics. Section VI presents conclusions and 
recommendations . 



3 







4 



IX. OVERVIEW OF TROPICAL PACIFIC OCEAN CIRCULATION: 

MEAN AND VARIABILITY 



To familiarize the reader with the area of the ocean 
studied, and as a basis for understanding later results, an 
overview of the observed mean circulation and its variability- 
in the tropical Pacific Ocean (20°S to 20°N) is presented. 

The mean zonal surface circulation of the tropical 
Pacific Ocean is characterized by alternating bands of 
eastward- and westward- flowing currents (Figure 1) which are 
forced by the large-scale wind field. The North Equatorial 
Current (NEC) and South Equatorial Current (SEC) form the 
equatorward arms of the subtropical gyres, and flow westward 
in the general area of the northern and southern Trades 
(Figure 2), respectively. Between these two currents (3°N to 
10°N) , the eastward flowing North Equatorial Countercurrent 
(NECC) is observed to flow counter to the prevailing winds. 
Also flowing eastward is the subsurface Equatorial 
Undercurrent (EUC) . The major part of the NEC curves to the 
northwest at about 150°E, past the Phillipines and Taiwan, to 
eventually feed the western boundary, the Kuroshio Current. 
The SEC, extending from about 3°N to 10°S, flows from South 
America to the western Pacific; as it approaches Australia, 
the SEC bifurcates near 18°S, with the southern branch 
feeding the East Australia Current and the northern branch 
continuing through the Solomon Sea, crossing into the 
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Northern Hemisphere to feed the NECC and the EUC in the 
western tropical Pacific (Figure 3) . A weak, eastward 
flowing South Equatorial Countercurrent (SECC) usually exists 
west of the International Dateline (IDL) at about 10°S. 

The NEC and SEC are directly wind-driven and are 
therefore strongly seasonal, reaching maximum velocity during 
winter in their respective Hemispheres. In the mid- Pacific, 
the major part of the NEC lies between about 8°N to 20°N, and 
carries about 45 Sv with a maximum speed of about 0.3 m/s in 
February. The SEC is strongest in August, at about 0.6 m/s, 
with a transport of about 27 Sv, decreasing to 7 Sv in 
February. The NECC varies seasonally in strength and 
position. From February to April the northwest monsoon 
prevents the SEC from feeding the NECC, so that the NECC is 
restricted to 4°N to 6°N with a volume transport of 15 Sv and 
maxiumum speeds below 0.2 m/s. East of 110°W the NECC 
disappears altogether. From May to January the NECC flows 
between 5°N and 10°N with speeds of 0.5 m/s. The annual mean 
transport of the NECC decreases uniformly with longitude, 
from 45 Sv west of 135°E to 10 Sv east of the Galapagos 
Islands . Maximum velocities of the SECC occur during the 
Northwest Monsoon (the cause of the wind stress minimum, 
February to April) . Typical surface speeds are below 0.3 m/s 
at 170°E, yielding a transport of about 10 Sv. In both 
Monsoon seasons the strength of the SECC decreases rapidly 
east of the IDL, and may be altogether absent from the 
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eastern Pacific (Tomczak and Godfrey, 1994) . 

Mean meridional circulation in the tropical Pacific is 
driven by poleward Ekman drift in the surface layers and 
equatorward motion at the depth of the thermocline. This 
circulation implies equatorial upwelling which can be 
measured from the divergence of the horizontal currents as 
measured with satellite tracked drifters (Philander, 1990) . 

With the exception of a band of relatively warmer 
surface water in the eastern tropical Pacific, just north of 
the equator, mean sea surface temperatures (SST) are higher 
in the western tropical Pacific than in the east. The 
eastward flow of the NECC coincides with this band of warmer 
surface water in the eastern Pacific.. SST in the eastern 
tropical Pacific are at a maximum during the northern spring, 
when the trade winds relax, and are at a minimum during the 
northern summer and fall when the trades are strongest (Apel, 
1987) . The expanse of very warm water in the western 
equatorial region has annual mean values exceeding 29.5°C 
between New Guinea and the Solomon Islands. Annual means 
above 29. 0°C. occur between about 10°N and 5°S, in the 
vicinity of the Solomon Sea, and extend eastward to about 
170°W. The thermal equator follows the Intertropical 
Convergence Zone (ITCZ) from 5°N in the west to about 15°N in 
the east, with annual means of about 27°C. (Tomczak and 
Godfrey, 1994) . 

Variability relevent to this study spans temporal scales 
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from intraseasonal to interannual. The interannual 
peturbation of the ocean-atmosphere known as the El Nino 
Southern Oscillation (ENSO) is a major component of the 
tropical Pacific variability, and is initiated by reversals 
of the easterly Trade winds in the western Pacific. These 
westerly wind bursts produce reversals of the westward 
flowing SEC across the basin. West of the dateline the 
reversals are in response to direct wind forcing while to the 
east they are caused by the passage of remotely wind-forced 
equatorial Kelvin waves with phase speeds of 2-3 m/s 
(McPhaden, 1993). These surface flows produce an eastward 
expansion of the western Pacific warm pool; the 
redistribution of warm water zonally flattens the east-west 
slope of the thermocline, weakening the baroclinic zonal 
pressure gradient and the EUC driven by it. ENSO tends to be 
phase locked to the annual cycle, with the largest SST 
anomalies occurring over the. equatorial cold tongue region in 
the month of December (Yang et al . , 1977). 

The intraseasonal variability addressed here spans 
temporal scales of 20 to 100 days, typified by barotropic and 
baroclinic instabilities associated with shear zones. The 
weakening or disappearance of the EUC impacts the generation 
of 20- to 3 0-day barotropic instability waves that are 
normally generated in the shear zone between the northern 
flank of the EUC and the SEC. These waves have westward 
phase velocity, wavelengths of about 1000 km, and have been 
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documented along 5°N in the shear zone between the SEC and 
the NECC by Perigaud (1990) using Geosat data. Westward 
propagating waves were documented with periods of 50-90 days 
and wavelengths of 630-950 km, along 12°N, in the shear zone 
between the NECC and the NEC, using the same data set. 
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Figure 1. Schematic of the surface currents of the Pacific 
Ocean. Abbreviations used in the area of study: Mindanao Eddy 
(ME) , Halmahera Eddy (HE) , and New Guinea Coastal Current 
(NGCC) . Figure 8.6 from Tomczak and Godfrey (1994) . 
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CURRENTS TRANSPORTS WINDS 




Figure 2 . Schematic of the relationships between the wind 
direction, surface current, and the Ekman transport (total 
transport in the wind-driven layers , shown by short arrows ) , 
in equatorial latitudes. Figure 5.1(a) from Ocean 
Circulation, G. Bearman, ed. , 1996. 



11 




120*E 140*E 160*E 



Figure 3 . Schematic of the surface currents of the western 
tropical Pacific. Abbreviations used: North Equatorial 
Current (NEC) , South Equatorial Current (SEC) , North 
Equatorial Countercurrent (NECC) , South Equatorial 
Countercurrent (SECC) , Halmahera Eddy (HE) , and Mindanao Eddy 
(ME) . Possible sources of the Indonesian Throughflow (IT) 
are shown as dotted lines . Figure 1 from Godfrey et al . 

(1993) . • 
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III. MODEL OUTPUT AND DRIFTER DATA SET 



A . MODEL DESCRIPTION 

Model fields are obtained from the LANL 1/6° global POP 
simulation, an eddy- resolving, z-level primitive equation 
model with active thermohaline dynamics (Maltrud et al., 

1997) . Smith et al. (1992) restructured the global, 1/2° 
Semtner and Chervin (1988) model derived from Bryan's (1969) 
formulation and versions ‘by Cox (1970, 1984) and Semtner 
(1974) , for a massively parallel Connection Machine 5 at LANL 
and configured it for a 1/6° grid. Dukowicz and Smith (1994) 
added an implicit free-surface formulation of the baro tropic 
mode. A full description of the model formulation can be 
found in Maltrud et al. (1997); here only relevant details 
will be given. 

1. Computational Grid 

To better track eddy length scales in this high 
resolution global simulation, a Mercator grid with 1280 
longitudinal and 896 latitudinal points was employed. The 
longitudinal spacing. was 0.28° and the latitudinal spacing 
equals the longitudinal spacing times the cosine of the 
latitude, producing a square grid over the global domain 
between 77°S to 77°N. The resolution ranged from about 31 km 
at the equator to just over 6.5 km at the highest latitudes, 
which on average, is 0.17 or 1/6°. The model has 20 non- 
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uniformly spaced vertical levels; the uppermost is 25 m deep, 
increasing to 550 m at the lowest level. 

2. Initialization and Forcing 

POP was initialized using output interpolated from the 
end of a 3-year Semtner (1997) 1/4° run (January 1988) , which 
in turn, had been initialized with results from an 
integration of 35 years from Levitus (1982) at 1/2° and a 
further 5 -year equilibrium run at 1/4°. POP then 
equilibrated for 5 years at the higher resolution, after 
which a 9-year run was completed. A density field from the 
end of this first POP run was used to initialize POP7 , the 
simulation used in this thesis. POP7 was forced with 3 -day 
averaged wind stresses derived from European Center for 
Medium Range Forecasting (ECMWF) winds from 1985-1994 (1000- 
mb winds from 1985-1989 and 10-m winds from 1990-1994) . 
Temperature and salinity in the upper level of the model were 
restored to Levitus (1982) climatology with a restoring 
timescale of 30 days. Very realistic topography (ET0P05) was 
employed. Biharmonic horizontal friction whose coefficients 
have a latitudinally decreasing form to take advantage of the 
square Mercator grid, and vertical mixing according to 
Pacanowski and Philander (1981) were used. 
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3 . Output Fields 

Extracted model fields are 3 - day snapshots of zonal 
velocity (u ) , meridional velocity (v) , and temperature (T) , 
for the period September 1992 to October 1994. High- 
frequency velocities were available for this period alone; 
only monthly averages were saved for the rest of the run. 
Fields were extracted from 112°E to 72°W, 21°S to 21°N, to 
coincide with the region from which the Pacific Ocean surface 
drifter measurements were collected. The resultant analysis 
grid has 655, 152, and 1 cells in longitude, latitude, and 
depth, respectively. Model output is filtered and binned 
into 2° latitude x 8° longitude areas for comparison with the 
drifter data. 

4 . Inertial Oscillations and Aliasing 

Forcing of high-resolution primitive-equation models 
with realistic, high frequency (period order of days) wind- 
stress fields produce inertial oscillations. Subsampling 
model output at intervals of several days leads to aliasing 
of these inertial currents into lower frequencies that vary 
with latitude (Jayne and Tokmakian, 1997) . Three-day wind- 
stress averaging and the use of 3 -day snapshots both 
contributed to the aliased signal in the POP output. 

Replacing the daily winds averaged over three days with daily 
winds interpolated to every timestep, and using three-daily 
averages rather than three-daily snapshots increased the 
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model energy levels (R. Tokmakian, personal communication) . 
Saving values of the prognostic variables at every time step 
and then averaging every three days avoids aliasing inertial 
oscillations (Jayne and Tokmakian, 1997) . Running another 
simulation, however, is not realistic, so the application of 
a 9 -day running mean is used to make the ouput serviceable, 
removing much of the spurious energy associated with the 
inertial motions. Figure 4 (R. Tokmakian, personal 
communication) shows the power spectrum of the u component of 
velocity, from the original and averaged runs, sampled at the 
local minima of the forcing function at 14 . 84°N and 14.84°S, 
along 160°W. At these latitudes, the period of the inertial 
oscillations is aliased to a period longer than the Nyquist 
frequency of the output data. Most noticeable is the 
difference in energy between the 3 -day snapshot and 3 -day 
average at periods less than 10 days. To remove much of this 
aliased signal, instantaneous 3-day snapshots were filtered 
using a 9 -day running average. Figure 5 (R. Tokmakian, 
personal communication) shows the reduction in amplitude of 
both signals at periods of less than 10 days, with much of 
the aliased material removed, after application of the 9-day 
running filter to the 3-day snapshots. 

5 . Previous POP Studies 

Previous relevant investigations and validation studies 
have shown the robustness of this model version. Fu and Smith 
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(1996) focus on the comparison of model results with TOPEX, 
while McClean et al . (1996) compare model results with TOPEX, 

as well as the Semtner-Chervin 1/4° model. Maltrud et al . 

(1997) examined mean circulation and variability in the 
model, with the effects of monthly varying and daily varying 
winds . 

Fu and Smith (1996) studied comparisons of the mean 
circulation, the mesoscale variability, the amplitude and 
phase of the annual cycle, as well as intraseasonal and 
interannual changes, and showed that the simulations and 
observations agree fairly well over a broad range of spatial 
and temporal scales. The sea level variance, however, 
produced by the model is generally less than the model by a 
factor of 2, primarily in the eddy-rich regions. The authors 
conclude that even higher resolution is needed to fully 
resolve the mesoscale eddies . 

Maltrud et al. (1997) state that many features of the 
wind-driven circulation are well simulated by POP and are 
relatively independent of various surface forcings applied. 
Surface current patterns from sea surface mean heights , and 
overall transports from the stream function of volume 
transport, are very realistic over most of the ocean, with 
the strongest simulated currents having width scales and core 
speeds comparable to their real counterparts. The authors 
cite very good agreement with observations of the Antarctic 
Circumpolar Current and the Pacific Equatorial Undercurrent, 
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and the data in the equatorial Pacific regime in general. Z- 
level models have typically had problems with the equatorial 
undercurrent due to vertical shear, but a meridional section 
of the zonal velocity at 150°W (Figure 6), shows the 
undercurrent and alternating jets decreasing in strength with 
depth. This good agreement in vertical shear, where previous 
simulations were problematic, as well as in horizontal 
circulation, are encouraging features in the model simulation 
of the equatorial Pacific. 

McClean et al . (1996), in descriptions of the spatial 

distribution of energy relative to T/P, concentrated upon 
statistical intercomparisons of mesoscale variability of T/P 
altimeter data, the 1/4° Semtner-Chervin model, and the 1/6° 
LoANL POP model; independent comparisons were made with 
drifter data in the tropics . Model results are in rough 
agreement with the T/P field in terms of spatial variability, 
with the best agreement in magnitudes in areas of strong 
currents, such as western boundary currents. Moderate 
variability is seen several degrees north and south of the 
equator in the Pacific, attributed to 20- to 30-day 
instability waves. Between 15°N and 15°S, drifter rms 
velocity field values are in close agreement with POP values, 
with drifter values exceeding total model velocities by a 
factor of 2 in the mid- latitudes . 
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B . DRIFTER DATA SET 

Under the sponsorship of the Tropical Ocean and Global 
Atmosphere (TOGA) Program of the World Climate Research 
Program (WCRP) , over 1300 satellite-tracked surface drifters 
were released between January 1979 and May 1994 in the 
tropical Pacific (20°N to 20°S) to observe surface currents 
and SST patterns across the entire basin (Niiler et al., 

1997) . Figure 7 indicates the frequency distribution of the 
number of five-day drifting buoy velocity observations as a 
function of year, showing the concentration in the years 
1992-1993. Later in this time period, these data include 
instruments deployed as part of the WOCE Surface Velocity 
Program. Today the WOCE and TOGA programs and their follow- 
on programs using similar drifters are referred to as the 
Global Drifter Program (GDP), managed by NOAA/AOML in Miami. 

1 . Satellite-Tracked Drifter Description 

The GDP drifters were drogued to 15 m depth in order to 
occupy the middle of the summer mixed layer and to lie below 
the daily thermocline and associated rotational currents, and 
were tracked to within a 300 m radius by Service ARGOS. A 
linear SST sensor accurate to 0 . 1°C is located on the lower 
half of the surface float. The drifter is fitted with a 
submergence sensor which allows the drogue status to be 
determined (drifter floats rarely submerge without drogues) . 
Drifter technical specifications and manufacturing procedures 
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are described in Sybrandy and Niiler (1991) . 



2 . Temporal and Spatial Extent of the Data 

The drifter data span the 1991-92-93 El Nino and the 
1988 La Nina; other variabilities include contributions from 
seasonal cycles, and eddy and equatorial wave motions of 
periods longer than 10 days. The data used in this study 
consists of 5-day-average drifter velocities binned to 2° 
latitude by 8° degrees longitude. 

Concentration of tracks across the tropical Pacific 
(Figure 8) is more a function of the circulation than 
location of release. Niiler et al . (1997) note that the 

measured circulation is the horizontal circulation of a 
three-dimensional flow, where the meridional and vertical 
velocity components play as significant a role as the mean 
zonal flows. Drifters deployed near the Equator generally 
moved to higher latitudes in response to mean Ekman 
divergences, which are not present in dynamic topography. 

Data density was not spatially uniform, as the drifters spent 
about twice as much time in convergent regions, such as the 
southern margin of the NECC, than at the Equator. 

3 . Sources of Error 

Error sources in constructing current fields from 
drifter data are from ARGOS position errors, the fidelity 
with which the buoys measure water movement (McNally, 1981) , 
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and the variability inherent in mesoscale perturbations 
(Hansen and Paul, 1984; Poulain and Niiler, 1989) . Niiler et 
al. (1995) established the drifter slip for tropical, GDP- 
style drifters, to be within 1 cm/sec in winds of 10 m/sec; 
drifter data used here have not been corrected for wind 
induced slip because the mean wind speeds in the tropical 
Pacific are less than 10 m/sec (Niiler et al . , 1997). Service 
ARGOS has consistently provided excellent positioning, though 
expensive, with uncertainty of less than 500 m (Poulain and 
Niiler (1989); this data set has an estimated accuracy of 300 
m (Niiler et al . 1997). The greater uncertainty associated 
with these data is that drifters do not provide a random 
sampling of currents since drifter deployment is not random 
in space or time, and the mean current near the equatorial 
divergence creates data sparse regions (Niiler et al . , 1997; 
Poulain 1993; Reverdin et al . , 1994). This data error is not 
estimated here, but we investigate this issue by comparing 
different sets of model-derived trajectories where we can 
control the initial distribution. 

Since water parcels in a non-dif fusive medium move on G_ 

surfaces and surface drifters on isobaric surfaces, 
trajectories of buoyant floats are not strictly Lagrangian 
(Fieux et al . , 1994; Freeland et al . , 1975; Poulain and 
Niiler, 1989) . Diffusivity effects can induce uncertainty, 
particularly in drifter deployments of long duration. 
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Uncertainties also occur in interpretation of mean flow 
information from ensemble averaged float velocities. Time- 
averaged flow over large areas is equal to the instantaneous 
average over many widely dispersed floats, if the turbulence 
is homogeneous and stationary. But if the eddy intensity 
varies in space, the center of mass of the float cluster will 
ascend the energy gradient to yield an inaccurate measure of 
mean flow (Freeland et al., 1975). This leads to apparent 
mean velocities if data concentration is not uniform and has 
been termed array bias (Davis et al . , 1996) . 

4 . Relevant Drifter Studies 

The use of surface floats for measuring current systems 
has become increasingly common as an effective method to 
measure surface circulation, as a complement to other 
observations, and as a method to refine numerical models with 
direct observations. A selection of drifter studies relevant 
to the area of study, model output domain, analysis methods, 
and dynamics involved is outlined here . 

Investigations into Pacific surface flow using 
Lagrangian surface floats include Hansen and Paul (1984) , 
with estimates of equatorial divergence in the eastern 
tropical Pacific from drifter-derived measurements of surface 
currents, and McNally (1981) in the northeastern Pacific. 
McNally et al . (1983) used the trajectories of 16 satellite- 

tracked drifting buoys and derived velocities to construct a 
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Lagrangian perspective of the near-surface circulation of the 
North Pacific subtropical gyre. Paduan and Niiler (1993), 
using Lagrangian surface drifters, report on the general 
nature of the large spatial scale and seasonal timescale 
currents and temperatures obtained from the drifters, and 
calculate the large-scale spatial convergence of heat in the 
northeast Pacific Ocean. Poulain (1993) presents estimates 
of horizontal divergence and vertical velocity in the 
equatorial Pacific, and Poulain and Niiler (1989) describe 
the statistical analysis of drifter motion in the eastern 
Pacific. Ralph and Niiler (1995) discuss preliminary results 
of the wind-driven currents deduced from Pacific drifter data 
sets, including the data set used here. Reverdin et al. 
(1985), Niiler et al. (1989), and Richardson (1983) use 
drifters to describe mesoscale eddies. Niiler et al . (1997) 

describe the data set used in this study. 

Hansen and Paul (1984) used data from 20 satellite- 
tracked drifting buoys deployed in the eastern equatorial 
Pacific in a study of the nature, effects, and energetics of 
currents associated with cusp-shaped long waves. They found 
that long waves are associated with a vigorous pattern of 
mesoscale eddies lying between the equator and 7°N, between 
110°W and 105°W. 

Poulain (1993) presents estimates of horizontal 
divergence and vertical velocity in the equatorial Pacific 
calculated from 12 years of surface drifter trajectories. 
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Mean divergence estimates are predominantly meridional and 
are similar to values from previous studies. Seasonal 
variability of equatorial divergence and upwelling are 
substantial and in phase with local wind forcing. 

Poulain and Niiler (1989) present a description of the 
surface circulation off the Southern California and Baja 
coasts, based upon statistical kinematic analysis of the 
trajectories of satellite tracked drifters. Lagrangian and 
Eulerian statistics of the variations about a mean southward 
drift are computed, with drifter velocities averaged in 200 
km by 200 km spatial bins. Absolute and relative particle 
dispersion was investigated using single and two-particle 
Lagrangian statistics. 

A preliminary report by Ralph and Niiler (1995) presents 
findings of wind-driven currents described from this set of 
satellite-tracked drifters in the tropical Pacific. Wind- 
driven currents were inferred by subtracting the geos trophic 
component from the observed total current; the residual or 
Ekman currents can be used to test models of vertical mixing 
and heat flux. Five-day segments of drifter velocity 
measurements were binned into 1 latitude degree x 5 longitude 
degree boxes and averaged within each box. The drifter data 
were used to assess the validity of the Ekman and Pollard 
models by comparing the magnitudes of the observed and the 
predicted velocities . 

Reverdin et al . (1994) use buoy drifts and current meter 
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records to present a statistical analysis of the seasonal 
variability of the equatorial Pacific Ocean at 15 meters. 

The objective of the analysis is to reduce the sampling noise 
by having a large data set, with the question of whether all 
drifters are comparable and whether the differences with the 
true currents can be corrected. The seasonal cycle is 
compared to ship-drift data, confirming that the analysis 
captured the seasonal cycle in the currents, and that the 
data is a faithful description of the mean tropical Pacific 
surface currents. The seasonal cycle is analyzed between 
20°N and 20°S on a 1° latitude by 5° longitude grid, which 
captures a large zonal seasonal variability of the currents 
within 15 degrees of the equator. 
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Figure 4. Power spectra (cpd) of the u velocity component at 
14.84°N and 14.84°S, along 160°W in the tropical Pacific, 
before application of the 9-day running filter. (R. 

Tokmakian, personal communication) 
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Figure 5 . Power spectra (cpd) of the u velocity component at 
14.84°N and 14.84°S, along 160°W in the tropical Pacific, 
after application of the 9-day running filter. (R. 

Tokmakian, personal communication) 
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Zonal Velocity at 150 W, Run 11 




Latitude 

Figure 6. Cross-section of the zonal velocity at 150°W from 
POP11 . The contour interval is 2 cm sec-i for u less than or 
equal to 10 cm sec* 1 , and 10 cm sec-1 for u greater than 10 cm 
sec-i. Dashed contours denote westward flow. Figure 10 from 
Maltrud et al., 1997. 
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Figure 7. Frequency distribution of the number of five day 
drifting buoy velocity observations as a function of year. 
(P. Niiler, personal communication). 
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Figure 8. Tracks of drogued drifters released in the tropical 
Pacific from January 1, 1988 to December 31, 1994. From 
Niiler et al. (1997) . 
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IV. 



EULERIAN RESULTS 



A. MEAN FIELDS 



1. Mean Zonal Velocity 

The mean zonal velocity (cm s -1 ) model and drifter fields 
(Figure 9) reveal the conspicuous surface circulation systems 
of the tropical Pacific, showing the low-latitude limbs of 
the recirculating subtropical gyres. Both data and output 
show the westward directed North Equatorial Current (NEC) , 
between 20°N and 8°N, the South Equatorial Current (SEC) , 
between 3°N and 10°S, and between them, the eastward flowing 
North Equatorial Countercurrent (NECC) . 

A weak eastward equatorial flow with values of 10 cm s -1 , 
west of the International Dateline and centered along the 
equator, is also common to both fields. Other similarities 
include a SEC maximum at 2°N, 125°W of about 50 cm s -1 ; NEC 
current values ranging from about -20 cm s -1 at approximately 
12°N and 140°E, to -10 cm s -1 in the eastern Pacific; and NECC 
maxima in excess of 30 cm s -1 in the western basin (5°N, 

130°E) , and in the eastern basin (5°N, 100°W) , with 
consistent values of at least 10 cm s" 1 across the basin. 

In addressing features expected from historical 
observations, Lukas (1987) also reports a SEC maximum of 50 
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cm s 1 at 2°N, and a minimum on the equator, along meridional 
transects at 158°W, 153°W, and 150°W. Drifter data would 
indicate a weak eastward flow along the equator at 135°W, 
previously observed by Hansen and Paul (1984) in 
investigations between 8°N and 8°S, 105°W to 110°W, and cited 
as imposing a jetlike structure in the westward flowing SEC. 
No indication of this weak eastward flow is indicated in the 
model field. As mentioned previously, the eastward flow in 
the western equatorial Pacific is indicative of westerly wind 
bursts . 

The most notable difference in the model field is the 
elongation of the SECC to the eastern basin, from 150°E to 
110°W. The latitudinal extent of the NEC in the model is 
slightly greater than that of the drifter data, extending 
from 5°N to slightly above 10°N, while the drifter data 
extends from 5°N to just below 10°N. (Reverdin et al., 1994). 

This discrepancy between model and drifter results may 
be due to the representation of the winds forcing the model 
in the tropical South Pacific; if there is a paucity of wind 
observations in this area, it is likely that the forcing 
could be somewhat erroneous . 

2 . Mean Meridional Velocity 

The most significant feature of both mean meridional 
velocity fields (cm s" 1 ) is the cross -equatorial symmetry of 



32 



generally northward flow north of the equator and southward 
flow south of the equator (Figure 10) . Maximum northward and 
southward flow for the drifter field is at approximately 
110°W, at 8 cm s _1 north of the equator and -10 cm s" 1 south 
of the equator. Model field maxima are displaced west to 
about 140°W, which may again be due to inadequacies in the 
wind forcing. 

Similar observations from drifter data were made by 
Luther and Johnson (1990) , and Hansen and Paul (1984) ; both 
studies report that, with the exception of brief crossings in 
conjunction with eddies, drifters avoided crossing the 
equator. Hansen and Paul (1984) infer from the divergence of 
the meridional flow with the equator that the dominant 
process is Ekman transport, divergent at the equator and 
driven by the southeast trade winds. The feature at 5°N, 

124°E (positive / negative reversal) in both model output and 
drifter data indicate a standing Rossby wave, where the mean 
flow has arrested the westward propagating Rossby wave and 
rendered it stationary (Lukas, 1987; Philander, 1990; 

Pedlosky, 1996) . In a study of the meridional structure of 
near-equatorial currents, Poulain (1993) concluded that the 
zonally averaged divergence is concentrated in a narrow 20-km 
band centered on the equator, and that mean divergence 
estimates are predominantly meridional. 

Significant differences exist in magnitudes, with the 
model output showing ranges in the central basin in excess of 
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18 cm s -1 at 4°N 140°W and -16 cm s -1 at 3°S 130°W, while the 
drifter data values are displaced to the east, with values in 
excess of 8 cm s -1 at 3°N 115°W and -10 cm s -1 at 3°S 115°W. 

These differences in magnitudes may be due to a lack of 
a model mixed layer. Since the model is solely dependent on 
the Richardson-number formulation of Pacanowski and Philander 
(1981) for the vertical redistribution of properties, the 
Ekman drift is largely trapped in the upper layer, producing 
artificially strong flow. 

3. Mean Velocity Vectors 

The velocity vector fields (cm s -1 ) again show the 
alternating eastward and westward zonal flows of the tropical 
Pacific (Figure 11) , with the latitudinal shear between 
oppositely directed currents. The southern branch of the SEC 
(centered at about 5°S) extends across the basin to New 
Guinea while the southern branch ends near the IDL. The East 
Australia Current is seen south of 15°S, along the 155°E 
meridian. 

In plots of surface current vectors of the tropical 
Pacific Philander et al. (1986) and Reverdin et al . (1994) 

also show the latitudinal shear north of the equator and 
strong westerly flow in the eastern Pacific. Reverdin et al . 
(1994) also show the bifurcation of the SEC in the western 
Pacific. 
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The most obvious difference between output and data 
occurred in the vicinity of the Indonesian Throughflow ( ITF) , 
with the highest drifter values indicating a strong cyclonic 
circulation in the Banda Sea that was not indicated in the 
model output; however, with the large grid size and few 
drifters in this region, no definite conclusions can be 
drawn. 



4 . Mean SST 

Mean temperature fields (°C) from model and drifter data 
exhibit the characteristic warm pool (designated by the 29°C 
isotherm in the drifter data) extending from the western 
equatorial Pacific southeast to about 160°W, 10°S, and the 
cold tongue extending along the equator from the South 
American coast to the mid-basin (Figure 12) . Except for a 
band of warm surface water that extends across the ocean 
basin just north of the equator, SST from both fields are 
much lower in the eastern than the western tropical Pacific. 
Hansen and Paul (1984) report similar SST patterns in the 
eastern equatorial Pacific and find evidence of cold 
intrusions that associate this cold tongue with equatorial 
divergence of Ekman transport and/or advection from the 
Peruvian coastal upwelling region to as far as 160°W. 

The drifter temperature field shows the western warm 
pool as larger and extending farther west than depicted by 
Levitus (1982) . Since the model is restored to Levitus, 
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conversely the western warm pool is not as extensive in the 
model output as it is in the drifter data. 



B . VARIABILITY 

The variability of the velocity components and 
temperature were calculated from the deviations of the bin- 
averaged ensemble means discussed in the previous section. 

1. Variability of Zonal Velocity 

Both fields show coherent patterns across the Pacific 
interior, with weaker variability poleward of 10°N and 10°S 
than in the equatorial region (Figure 13) . Drifter and model 
both show two general regions of relative maxima along the 
equator and at 7°N, and two general regions of relative 
minima, at 14°N and south of 10°S, east of 140°W. Model 
maxima are about 15 cm s -1 at 7°N and 25 cm s'" 1 at the 
equator, while drifter maxima for these regions are 25 and 50 
cm s -1 , respectively. Both fields have similar locations of 
maxima in the the western Pacific, at about 130°E and 4°N, 
with values of about 40 cm s _1 for drifter data and 30 cm s -1 
for the model values . 

Niiler et al . (1997) state that the maximum in the 

eastern basin on the equator, between 100°W and 140^, 
results from variations produced by equatorial wave 
propagation from the west, and from the local generation of 
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instability waves in the east. A second relative maximum is 
at about 7°N in the eastern basin, where organized wave/eddy 
motions have been reported in Geosat sea surface height data 
(Perigaud, 1990) , or possibly due to seasonal variability of 
the NEC. Poulain (1993) also mentions that the variability 
not associated with the wind could be induced by inertia- 
gravity waves excited by shear instability of the zonal 
currents. Reverdin et al . (1994) report high frequency (20 

to 30 days) values on the equator ranging from 14 cm s _1 at 
170°W to 26 cm s" 1 at 110°W. 

Model values are less than drifter data in all regions; 
along the equator the model western basin maximum displays 
about 80% of the drifter variability and is displaced to the 
west (possibly as a result of the wind forcing) . The 
magnitude of the eastern basin model maximum is about 55% of 
the observed value, while the model maximum at 7°N is about 
60% of the corresponding observed value. Away from these 
maxima, the model values are roughly half the drifter values. 

2 . Variability of Meridional Velocity 

Variability in v (cm s -1 ) for both fields is largest just 
north of the equator in the western and east-central Pacific, 
and north of the equator in the eastern Pacific (Figure 14) . 

In the western Pacific, drifter and model maxima are centered 
3°N, 130°E, with drifter values of >35 cm s -1 , and model 
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values of >14 cm s -1 . Maxima exist in both fields on the 
equator near 4°N, 128°W, with drifter values of >30 cm s" 1 , 
and model values of >10 cm s -1 . 

Reverdin et al . (1994), in using data from five current 

meter moorings from 1987 - 1992, report high frequency (20 to 
30 days) meridional variance values along the equator ranging 
between 15 cm s -1 at 170°W to 28 cm s -1 at 110°W. These 
values in the eastern Pacific are again caused by instability 
waves (Hansen and Paul, 1984; Reverdin et al., 1994). 

Poulain (1993) associates the maxima in the eastern Pacific 
with semiannual fluctuations of the equatorial divergence. 

The most significant difference exists in the 
magnitudes, with model output less than drifter values. This 
may be due to inadequate representation of wind forcing in 
the model, lack of a mixed layer, and possibly still too 
coarse resolution. 

3. Variability of Total (u and v) Velocity 

Both model- and drifter-based representations of eddy 
kinetic energy show good agreement (Figure 15) in the central 
basin straddling the equator, away from land, with drifter 
values (cm s" 1 ) slightly higher than model output; output and 
data do not agree as well at higher latitudes . Drifter data 
show greater variability in the vicinity of the HE and ME, in 
the south Pacific, east of Australia, and in the WBC region 
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north of New Guinea. Drifter variability is also greater in 
the region between 10°S and 20°S. 

4. Variability of SST 

Drifter data (Figure 16) show a temperature variance 
maxima (°C) along the equator, in the eastern Pacific 
extending from the South American coast to about 130°W, 
approximately at the IDL; a weaker maxima is seen in the 
model. Hansen and Paul (1984) also refer to a region of 
larger variance at 2°S, not indicated in the drifter or model 
data set, and attribute it to the passage of a cold pool, 
perhaps due to upwelling. Drifter variability is greater 
than model output due to model restoration to Levitus 
climatology, on a 30-day time scale. Also, applying a 9-day 
running mean to temperature, for consistency with the 
velocity data, further reduced the temperature variability. 



C . DIVERGENCE AND RELATIVE VORTICITY 
1 . Divergence 

Divergence (10" 7 s _1 ) was calculated using 
V»V W = (<?<«> I dx + d^v^ldy ) 

in a centered-dif f erence scheme applied to the 2° x 8° 
Eulerian means discussed above . 

Divergence in the drifter and model (Figure 17) was well 
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defined and due mainly to d<v>/dy. Both fields are similar in 
outline, with strongly divergent circulation along the 
equator east of the IDL and weakly divergent flow in the 
western basin. Both representations show equatorial 
divergence at all meridians, except at about 120°E. The 
highest value for the model output is 7 x 10~ 7 s' 1 , centered 
at about 140°W, while the drifter data has two regions of 
maxima, 4 x 10“ 7 s' 1 , centered at about 145°W, and 3 x 10" 7 s' 1 
centered at about 110°W. Both data and output indicate the 
southern limb of the convergent north Pacific subtropical 
gyre north of 12°N. The convergence, poleward of 3°N and 3°S 
in both model output and drifter data, was more broad scale 
than the divergent fields. The divergence is due principally 
to the surface poleward Ekman flow on each side of the 
equator, driven by the prevailing westward trade winds . 

The most significant discrepancy is slightly higher 
model magnitudes, due to use of mean velocity values in the 
divergence calculation. The divergence calculation depends 
mainly on d<v>/dy, which reflects the mean meridional velocity 
<v > . Thus, the generally higher mean velocity values from 
the model output, compared to drifter data, also translates 
to higher divergence values . 
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2. Relative Vorticity 

Relative vorticity (10~ 7 s" 1 ) was also calculated from 
the 2° x 8° Eulerian averages using 

C= (5<v>/A -d(uyijy), 

and is due principally to d<u>/dy, reflecting variations of 
the zonal currents. The relative vorticity of both model-and 
drifter-based fields (Figure 18) is similarly defined in 
spatial extent, with a narrow band of positive vorticity 
along the equator with negative vorticity north and south, 
alternating approximately every five degrees. Relative 
vorticity values are similar along the equator, with the 
model about 30% higher than observed data in the regions of 
10°S and .10°N. 

D . COVARIANCE FIELDS 

1. Covariance of u and v: <u'v'> 

Given N pairs of observations the covariance between two 
variables is calculated from (Chatfield, 1996) : 

s xy = E ( x i - *)()>« - 7 ) / (N - 1). 

i= 1 

Both drifter data and model output (Figure 19) show 
consistent negative values north of the equator. The drifter 
data show maxima along the equator, between the dateline and 
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100°W, associated with the velocity variance peaks. The high 
negative values in the western tropical Pacific of the 
drifter data occur farther to the west and north of the 
equator in the model; however, the drifter data is noisy. 

The covariance maximum along the equator in the drifter data, 
between the IDL and 100°W, is not indicated in the model 
output . 

Hansen and Paul (1984) indicate that the negative 
covariance north of the equator supports the hypothesis that 
waves derive energy from the barotropic instability of the 
meridional shear of the zonal mean flow. Luther and Johnson 
(1990) attribute the broad maximum along the equator to a 
momentum flux that is down gradient of the mean zonal 
velocity, with magnitudes comparable to a WBC. 



2. Covariance of Velocity and SST 

The covariance of zonal velocity and SST (cm s -1 °C) 
in drifter and model (Figure 20) vector plots both show a 
tendency of the covariance vectors to tilt towards the 
equator, implying a net eddy heat convergence in the 
equatorial band. Both exhibit rough similarity in the 
general distribution of the positive covariance along the 
equator, indicating eastward heat flux, with the model output 
amplitudes less than the drifter data. Both fields indicate 
convergence (downwelling) in cooler regions (along 5°N and 
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5°S) , where poleward flow predominates, and divergence 
(upwelling) in warmer regions (along the equator) where 
poleward flow is weaker. The greatest variance in both 
fields is shown in the cold tongue along the east basin, 
extending westward into the high gradient region of the 
tongue. The prevalence of equatorial upwelling is shown as 
this narrow tongue of relatively cooler SST extending from 
the coast of South America to the IDL. The drifter data 
values are highest in the east, at approximately 100°W, 
whereas the model output is greatest at the IDL with smaller 
but consistent values to the east. 

Hansen and Paul (1984) report a slight negative 
covariance in the southern hemisphere (as indicated in the 
drifter data, but not apparent in the model output), and more 
definite negative values north of the equator (indicated in 
the model but less apparent in the drifter data) . 
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Figure 9 . Mean Zonal Velocity (cm s _1 ) : Drifter Data (top) and 
Model Output (bottom) . 
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Figure 10. Mean Meridional Velocity (cm s-i) 
(top) and Model Output (bottom) . 
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Figure 11. Mean (u and v ) Velocity Vectors (cm s _1 ) : Drifter 
Data (top) and Model Output (bottom) . 
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Figure 12. Mean Temperature (°C) : Drifter Data (top) and 
Model Output (bottom) . 
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Figure 13. Variability (standard deviation) of Zonal Velocity 
(cm s-i) : Drifter Data (top) and Model Output (bottom) . 
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Figure 14. Variability (standard deviation) of Meridional 
Velocity (cm s-i) : Drifter Data (top) and Model Output 
(bottom) . 
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Figure 15. Variability (square root of eddy kinetic energy) 
of u and v Velocity (cm s-i) : Drifter Data (top) and Model 
Output (bottom) . 
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Figure 16. Variability of SST (°C) : Drifter Data (top) and 
Model Output (bottom) . 



51 



LATITUDE LATITUDE 




120”E 160*E 160°W 120°W 80°W 

LONGITUDE 



Figure 17. Divergence ( 10-7 s -i) : Drifter Data (top) and Model 
Output (bottom) . 
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Figure 18. Relative Vorticity (10 -7 s _1 ) : Drifter Data (top) 
and Model Output (bottom) . 
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Figure 19. Covariance of Total (u and v ) Velocity (cm2 s - 2 ) : 
Drifter Data (top) and Model Output (bottom) . 



54 



LATITUDE LATITUDE 




120'E 160°E . 1 60°w 1 20"W 80”w 




1 20°E 1 60°E I60°W 1 20°W 80”w 

LONGITUDE 



21 -° 

Figure 20. Covariance of Velocity and SST (cm s~i °C) : Drifter 
Data (top) and Model Output (bottom) . 
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V. LAGRANGIAN STATISTICS 



The Eulerian statistics presented previously do not take 
full advantage of the particle-path (trajectory) nature of 
the drifter data. Indeed, Lagrangian statistics that 
describe the decorrelation time and length scales, along with 
the strength of the turbulent dispersion (i.e., the 
diffusivity) , by following the horizontal particle paths has 
been a primary focus of open-ocean drifter experiments 
(Poulain and Niiler, 1989; Paduan and Niiler, 1993). 

The large spatial and temporal extent of the drifter 
data set makes it possible to examine Lagrangian statistics 
of the tropical Pacific surface circulation. Poulain 
(personal communication, 1997) has made preliminary 
calculations of these quantities for the drifter data. This 
section presents the Lagrangian autocorrelation, time and 
length scales, and diffusivity based on particle paths 
created from the model velocity fields, for comparison with 
the drifter-based results. Statistical methods used here 
have been presented by Davis et al . (1996), Poulain and 

Niiler (1989), and Paduan and Niiler (1993), among others. 

Model velocity fields were used to produce simulated 
Lagrangian trajectories for comparison with actual drifter 
trajectories. An advantage of model-based trajectories is 
the ability to initialize the trajectories at any time/space 
point. In this study, two sample initialization grids were 
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chosen. One released particles at 460 uniformly spaced 
locations in the model grid using the 2° by 8° spacing of the 
Eulerian comparisons described earlier. The second data set 
consisted of 280 trajectories initiated at the actual 
deployment locations of the WOCE/TOGA drifters during the 
period September 1992 to October 1994. In this non-uniform 
case, all trajectories were initialized at the start of the 
model period, as was also the case for the uniform grid. 
Trajectories themselves (Figure 21) were obtained by using a 
fourth-order Runge-Kutta method (Appendix) as used by Hofmann 
et al . (1991). Lagrangian autocovariance and diffusivity 

functions from individual trajectories were ensemble-averaged 
according to mean positions (Figure 22) in latitude bands of 
22.5°N to 22.5°S, in 5° bands. The resulting Lagrangian 
statistics are described below and are listed in Table I and 
Table II. 

In a comparison of trajectories of the Pan-Pacific 
Surface Current Study (PPCS; Figure 8) and the uniform and 
non-uniform data (Figure 21) , the most significant features 
in common are the concentration of tracks in the western 
Pacific, at around 140°E, between the equator and 5°N in the 
uniform and non-uniform fields, and between the equator and 
5°S in the PPCS. Also, both model- and drifter-derived 
trajectories exhibit a paucity of tracks in the central basin 
along the equator. Track density is a function of deployment 
strategy and mean circulation patterns, primarily equatorial 
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divergence. This is made clear by the comparison of the 
uniform and non-uniform data sets in Figure 21. Although the 
very high concentration of trajectories in the western 
Pacific is largely the result of large numbers of deployments 
in that region, the uniformly deployed data set shows there 
to be convergence in the model velocities in that region. 
Similarly, the additional equatorial deployments in the non- 
uni form data set are not enough to overcome the strong 
divergence in that region. 



A. 



LAGRANGIAN AUTOCOVARIANCE 



The Lagrangian autocovariance is defined as the average 
time-lagged covariance of velocity along a trajectory (Davis, 
1993) : 



1 



r o +r 



R jj(x,T)X q, t^) - rj, J u -(t, x q, u j(t~¥ t, Xq, t^) dt i u,v,j u, v 



where u' is the perturbation velocity and T is the length of 
the time series. The velocity components u and v refer to 
the velocity at time t of the drifter passing through x Q at 

the initial time t 0< If the velocity field is assumed to be 

stationary and homogeneous, as was done in this study, then 
dependence upon the initial point x Q vanishes and the equation 
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simplifies to: 

1 7 

Rij(% T) = y $ u i'(t) Uj'(t+T)dt i = u,v; j = u,v . 

The zero-lag autocovariances R (0,T) and R (0,T) are 

uu w 

the velocity variances in the east-west and north-south 
directions, respectively, for the time series of length T. 

Lagged autocovariance functions were computed for each 
trajectory and regional (zonal) statistics were estimated by 
ensemble averaging all functions whose mean position fell 
within 5° latitude bands. For ensemble averages, the 
individual autocorrelation functions were calculated by 
normalizing the autocovariance by the variance before 
ensemble averaging . 

Ensemble averages, by latitude band, of the 
autocorrelation functions are shown in the center panel of 
Figures 22 - 30 for the uniform distribution and Figures 31 - 
39 for the non-uniform distribution. The standard deviation 
of the mean autocorrelation function is shown for each time 
lag. For both the uniform and non-uniform ensemble-averaged 
time-lagged autocorrelation functions, v decays faster in all 
cases, but error bars indicate that the averages are not well 
defined. R should smoothly asymptote to zero in the presence 
of stationary, homogeneous turbulence; the degree to which 
the functions fail to asymptotically approach zero provides a 
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measure of the non-homogeneity of the field. Most likely, 
deviations from the expected functional forms reflects 
inadequate description of the mean currents and their 
horizontal shear (Bauer et al . , 1997) so that a single 
turbulent (mesoscale) time scale is not appropriate. 

The ensemble -averaged time-lagged autocorrelation 
functions based on the uniform trajectories are well-behaved, 
with some indications of multiple time scales in the v 
component of 7.5°S to 2 . 5°S and the 2 . 5°N to 7.5°N bins. The 
bin spanning the equator, from 2.5°S to 2.5°N, shows a more 
rapid decay of the v component than other bins . The 
ensemble-averaged time-lagged autocorrelation functions based 
on the non-uniform trajectories show similar behavior in 
similar bins. 

B. INTEGRAL TIME AND LENGTH SCALES 

The Lagrangian integral time and space scales are the 
time and distance over which a drifter remembers its path, 
and are defined by (Poulain and Niiler, 1989) : 




1 
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The zonal and meridional autocorrelation functions 



computed here asymptotically approach zero, but in some cases 
do not cross zero. Therefore, the choice of e-folding scale 
rather than zero crossings was made to provide consistent 
estimates of time and length scales. Furthermore, the method 
used to determine the time scale selected the closest lag to 
the integral e- folding time resulting in a discretization of 
the estimates into the 3 -day time steps of the trajectory 
data. 

The ensemble -averaged time and length scales for each 
latitude band, by uniform and non-uniform deployment, are 
shown in Figures 41 and 42, respectively, and listed in Table 
I and Table II. In the tables, the range estimates derive 
from the e- folding times of the ensemble averaged functions 
plus and minus one standard deviation. 

Zonal integral time scales for the uniform field vary 
between 18 days, for the southernmost and northernmost bands, 
and 51 days for the 12.5°N to 7 . 5°N band. Meridional time 
scales, conversely, vary between 6 or 9 days in the seven 
innermost bands, to 12 days at the northern- and southern- 
most bands. Poulain (personal communication, 1997), reports 
a zonal integral time scale (Figure 43) of about 10 days, and 
a meridional time scale (Figure 44) of about 5 days, 
throughout four seasonal, ensembles. Zonal and meridional 
integral time scales for the non-uniform fields reflect the 
same distribution as the uniform field, despite the large 
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differences in trajectories per band for several latitude 
bands . 

Average length scales of the uniform field for east-west 
and north-south drifter motions are 130 km and 29 km, 
respectively, and are consistently larger in the east-west 
than in the north-south direction. Length scales seem more 
sensitive to drifter concentration than time scales. Zonal 
length scales for the uniform field vary between 24 km in the 
southernmost band, consistently increasing to the north, to a 
high of 283 km in the 12 . 5°N to 7 . 5°N band. Meridional 
length scales vary between 14 to 20 km in the four 
southernmost bands, then between 28 and 35 km in the next 
four bands to the north, and reaches a maximum of 67 km at 
the northernmost band. 

C. DIFFUSIVITY 

The diffusivity K, which is defined as the time rate of 
change of the displacement covariance, is related to the 
integral of the Lagrangian autocovariance, in the stationary, 
homogeneous regime by (Poulain and Niiler, 1989) : 

1 1 

K u (t)= T J R u (r) dr. 
z o 

The diffusivities K and K are computed from the individual 

uu w 
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autocovariance functions without normalization, and then 
ensemble -averaged by latitude bands for each lag as for R 
(bottom panel, Figures 23 - 40; note that large variations in 
K necessitate different ordinate scales for some figures) . 

To represent the theoretical asymptotic values for K uu 

and K , values were selected at the integral times computed 

from the normalized R and R above. These diffusivity 

values are an approximation to the asymptotic solution at the 
time lag corresponding to the integral time scale. These are 
plotted versus latitude in Figure 45 and listed in Table I 
and II. Range values derive from the dif fusivities plus or 
minus one standard deviation for the time lags corresponding 
to the integral time scales . 

In the uniform data the zonal component of diffusivity 
K uu dominates the results with maximum values ranging from 0.2 

x 10 3 m 2 s -1 at high latitudes to 12 x 10 3 m 2 s _1 near the 
equator. Typical magnitude of the meridional diffusivity K 

ranges from about 0.1 x 10 3 m 2 s _1 at 20°S to 2.2 x 10 3 m 2 s _1 at 
20°N. Poulain (personal communication, 1997) report maximum 
values of the zonal component of diffusivity (Figure 46) of 
50 x 10 3 m 2 s -1 at high latitudes to 150-200 x 10 3 m 2 s -1 near the 
equator, and a typical magnitude of meridional diffusivity 
(Figure 47) of 10-20 x 10 3 m 2 s -1 . Bauer et al. (1997) report 
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K and K values in the SEC of 7 x 10 3 m 2 s _1 and 3 x 10 3 m 2 s _1 , 

uu w 

respectively; they report K and K ' values in the NECC of 15 
x 10 3 m 2 s -1 and of 4 x 10 3 m 2 s -1 , respectively. 

D. POLARIZATION 



A Lagrangian measure of the polarization or mean angular 
momentum J is given by the integral of the covariance 
functions (Poulain and Niiler, 1989) : 

t 

J(t)=f RJt)-RJt) dr 

o 

The direction of rotation of the trajectories can be 
determined from combinations of the covariance functions 
between zonal and meridional currents, which provides a 
characterization of the eddy field. Positive J indicates 
cyclonic (counterclockwise) rotation, while negative J 
indicates anticyclonic rotation. Values of J are shown in 
the top panel of the uniform (Figures 23 - 31) and non- 
uniform (Figures 32 - 40) averaged ensembles. Because these 
functions depend on the Lagrangian autocovariance functions 
and the assumptions of homogeneity and stationarity, values 
for times greater than the integral time scales should be 
ignored. For the long model trajectories used here and the 
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rather large areas that these statistics represent, there is 
no evidence for significant bias toward clockwise or 
counterclockwise rotations in the trajectories. 
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Figure 21. Resultant trajectories induced by POP model 
velocity fields from 460 uniformly spaced initial positions 
(top) and 280 non-uniformly spaced initial positions 
(bottom). Crosses are .initial positions. 
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Figure 22. Mean and initial positions for 460 uniformly 
spaced initial positions (top) and 280 non-uni formly _ spaced 
initial positions (bottom) . Crosses are initial positions, 
circles are mean positions. 
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Figure 23. Ensemble Averaged Uniform Deployments. Latitude 
Band 22 . 5°S - 17.5°S. Polarization (10 7 m 2 s-i) or measure of 
angular momentum J (top) ; autocovariance R uu (solid line, 
error envelopes) and i! w (dashed line, error bars) (center) ; 
diffusivity K uu (m 2 s -1 ) (solid line, error envelopes) and K-w 
(dashed line, error bars) (bottom) . 
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Figure 24. Uniform Deployments. Latitude Band 17.5°S - 
12 . 5°S . 
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Figure 25. Uniform Deployments. Latitude Band 12.5°S - 
7 . 5°S . 



,5 Uniform Bin 4: -7.5 to -2.5 






-4000 



0 10 20 30 40 50 60 70 80 

time lag (days) 



Figure 26. Uniform Deployments. Latitude Band 7.5°S - 
2 . 5°S . 
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Figure 27. Uniform Deployments. Latitude Band 2.5°S - 
2 . 5 °N . 
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Figure 28. Uniform Deployments. Latitude Band 2.5°N - 
7.5°N. 
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Uniform Bin 7: 7.5 to 12.5 






Figure 29 . Uniform Deployments . Latitude Band 7 . 5°N - 
12 . 5°N. 
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Figure 30. Uniform Deployments. Latitude Band 12.5°N - 
17 . 5°N. 
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Uniform Bin 9: 1 7.5 to 22.5 






Figure 31. Uniform Deployments. Latitude Band 17.5°N - 
2 2 . 5 °N . 
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Figure 32. Ensemble Averaged Non-Uniform Deployments. 

Latitude Band 22.5°S - 17.5°S. Polarization (10 7 m2 s-i) or 
measure of angular momentum J (top) ; autocovariance R uu (solid 
line, error envelopes) and R-w (dashed line, error bars) 
(center) ; diffusivity K uu (m2 s-i) (solid line, error 
envelopes) and K-w (dashed line, error bars) (bottom). 
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Figure 33. Non-uniform Deployments. Latitude Band 17 . 5°S - 
12 . 5°S . 
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Figure 34. Non-uniform Deployments. Latitude Band 12.5°S - 
7 . 5°S . 
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Non-Uniform Bin 4: -7.5 to -2.5 






Figure 35. Non-uniform Deployments. Latitude Band 7 . 5°S - 
2 . 5°S . 
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Figure 36. Non-uniform Deployments. Latitude Band 2 . 5°S - 
2 . 5°N. 
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Figure 37. Non-uniform Deployments. Latitude Band 2.5°N - 
7 . 5°N. 
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Figure 38. Non-uniform Deployments. Latitude Band 7 . 5°N - 
12 . 5°N . 
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Non-Uniform Bin 8: 12.5 to 17.5 






Figure 39. Non-uniform Deployments. Latitude Band 12 . 5°N - 
17 . 5°N . 
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Figure 40. Non-uniform Deployments. Latitude Band 17.5°N - 
2 2 . 5 °N . 
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Figure 41. Meridional profiles of zonally averaged uniform 
(dashed line) and non-uniform (solid line) time scales. 
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Figure 42. Meridional profiles of zonally averaged uniform 
(dashed line) and non-uniform (solid line) length scales. 



88 



O) 

CD 

t 

CD 

r^ 

CD 

CO 

DC 

LU 

\— 

U. 

DC 

Q 

Q 

LU 

z> 

o 

o 





to o to o up o in 



• i 






hiuon aarunvi 



Figure 43. Zonal integral time scale T uu (days) from drifter 
trajectories as a function of latitude (90°W to 160°W at a 
resolution of 2° latitude) for the four seasons of the year. 
(P. Poulain, personal communication) . 
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Figure 44. Meridional integral time scale T w (days) from 
drifter trajectories as a function of latitude (90°W to 160°W 
at a resolution of 2° latitude) for the four seasons of the 
year. (P. Poulain, personal communication) . 
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Figure 45. Meridional profiles of zonally averaged uniform 
(dashed line) and non-uniform (solid line) diffusivity 
(m2 s _1 ) , and number of drifters per zonal band. 
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Figure 46. Zonal diffusivity K uu (103 m 2 s -i) from drifter 
trajectories as a function of latitude (90°W to 160°W at a 
resolution of 2° latitude) for the four seasons of the year 
(P. Poulain, personal communication) . 
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Figure 47. Meridional diffusivity ( 1 0 4 m 2 s-1) from 
drifter trajectories as a function of latitude (90°W to 160°W 
at a resolution of 2° latitude) for the four seasons of the 
year. (P. Poulain, personal communication). 
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36 


15 
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6 


44.6 


23 .0 


to 
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18 


9 


89.2 
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42 


15 
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57.5 
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to 


11766 


618 


51 


9 


283.8 


29 . 9 


7 . 5°N 


16395 


676 


87 


12 


484.1 
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16.9 


to 


12428 


1356 


33 
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3 


42.2 
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6 
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45 


6 
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9 
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45.3 
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24 
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51 
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3 


18.7 
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18 


6 


56.1 
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60 


15 


186.9 


36.0 
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355 


67 


9 


3 


19.0 


5.1 


to 
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145 


18 


9 


38.0 


15.2 


17 ,5°S 


1367 


172 


57 


12 


120.3 


20.3 


17 . 5°S 


92 


82 


6 


6 


8.1 


8 . 1 


to 


219 


134 


18 


12 


24.2 


16.2 


22 . 5°S 


371 


188 


45 


21 


• 60.4 


28.4 



Table I. Diffusivity, time and length scales for the uniform 
drifter distribution by latitudinal bands. Bold values are 
means, lower and upper values are plus and minus one standard 
deviation, respectively. 



94 



! NON-UNIFORM 
iModel Field 


Diffusivity 


Time 


Scale 

| 


Length 


Scale 


Latitude 


Kll 


K22 


U 


V 


U 


V 


Band 


(m 2 s-i) 


(m 2 s- 1 ) 


(days) 


(days) 


( km ) 


(km) 


22 . 5°N 


539 


351 


6 


6 


18.7 


15.8 


to 


1610 


456 


24 


9 


74.6 


23.7 


17 . 5°N 


2354 


597 


42 


21 


130.6 


55.2 


17 . 5°N 


1621 


368 


15 


6 


53.4 


16.7 


to 


2899 


604 


33 


12 


117.6 


33.4 


12 . 5°N 


4868 


1111 


75 


45 


267.2 


125.2 


12 . 5°N 


15734 


. 1007 


24 . 


6 


213.3 


29.0 


to 


28869 


1007 


60 


6 


533.3 


29 . 0 


7 . 5°N 


33240 


1151 


87 


9 


773.3 


43 . 6 


7 . 5°N 


7661 


1151 


15 


6 


117.5 


30.6 


to 


13805 


1416 


33 


9 


258.6 


45 . 8 


2 . 5°N 


20630 


1635 


66 


15 


517.1 


76.4 


2 . 5°N 


6693 


512 


9 


3 


83.0 


14 . 2 


to 


12073 


634 


21 


6 


193.8 


28.3 


2 . 5°S 


15905 


634 


36 


6 


332.2 


28.3 


2 . 5°S 


1561 


350 


9 


6 


41.2 


19.4 


to 


2545 


350 


18 


6 


82.5 


19.4 


7 . 5°S 


4162 


364 


39 


9 


178.7 


29 . 1 


7 . 5°S 


2178 


774 


6 


6 


37.9 


25.1 


to 


5499 


774 


21 


6 


132.8 


25.1 


12 . 5°S 


7987 


2374 


48 


30 


303.5 


125.6 


12 .5°S 


1962 


367 


3 


3 


23.8 


11.1 


to 


7698 


757 


15 


9 


119.1 


33 . 3 


17 . 5°S 


•11767 


1037 


30 


15 


238.2 


55.5 


17 ,5°S 


361 


36 


18 


3 


28.8 


3.5 


to 


414 


126 


21 


18 


33.6 


21.2 


22 . 5°S 


797 


211 


54 


36 


86.4 


42.4 



Table II. Diffusivity, time and length scales for the non- 
uniform drifter distribution by latitudinal bands. Bold 
values are means, lower and upper values are plus and minus 
one standard deviation, respectively. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 



The results of the present study have underscored the 
complementary use of drifter observations and numerical 
simulations for studying ocean circulation. Model successes 
can be attributed to very realistic physics and forcing at 
high resolution. The model simulates the mean surface 
circulation of the tropical Pacific quite well, as evidenced 
by the mean Eulerian velocity fields, reproducing observed 
flow patterns. The discrepancies are typically in 
magnitudes, primarily in the equatorial regions. 

The model mean values exceed the observed values by as 
much as 50% in the equatorial region (Figures 9 and 10) . The 
higher model mean values occur probably as a result of a lack 
of a mixed layer. Since the model is solely dependent on the 
Richardson-nuxnber formulation of Pacanowski and Philander 
(1981) for the vertical transfer of momentum and tracers, the 
Ekman drift is largely trapped in the upper layer, producing 
an artificially strong flow. The unrealistic representation 
of the SECC may be due to an inadequate representation of 
wind forcing. 

Model variability (Figures 13 - 16) is about 20% less 
than the observed data in equatorial regions, and 50% less 
than the observed data poleward of 10°S and 10°N. The model 



97 



is not able to completely simulate variability in part due to 
a lack of an appropriate parameterization of mixed layer 
physics. The westward displacement of zonal variability 
maxima (Figure 13) may be due to inadequate representation of 
ECMWF winds in the western Pacific. Insufficient 
representation of variability has also been attributed to 
deficiencies in simulating the mean structure of the density 
field (Beckmann et al . , 1994), which can inhibit the 
generation of barotropic / baroclinic instabilities, 
associated with the horizontal and vertical shear of zonal 
currents. Applying the 9 -day running mean to the temperature 
field, for consistency with the velocity data, further 
reduced the temperature variability (Figure 16) . 

Relative vorticity (Figure 18) and divergence (Figure 
19) are dominated by d<u>/dy and d<v>/dy, respectively; 
therefore, the structure of these fields reflects the 
distribution of the respective dominant mean flow component. 
The generally higher mean velocity values from the model 
output, compared to drifter data, also translates to higher 
divergence and relative vorticity values, relative to 
observations . 

Both model- and drifter-based fields of u and v 
covariance (Figure 19) show generally negative values north 
of the equator and positive values south of the equator. The 
strong gradients in the eastern equatorial region, 
particularly in the north, are attributed to instability 
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waves. The covariance of velocity and SST (Figure 20) 
indicates convergence of heat toward the equator, with the 
strongest values just north of the equator, most likely from 
20-30 day instability waves. The model extends the strong 
equatorial fluxes farther west and more symmetrically about 
the equator. 

In comparing the Lagrangian results of the uniform 
distribution of model trajectories to observed values, model - 
based time scales are too slow by as much as a factor of 2, 
the length scales are as much as four times higher than 
observed values in the equatorial region (Reverdin et al., 
1994; Hansen and Paul, 1984), and the dif fusivities are too 
low by an order of magnitude. These differences can be 
attributed in part to the lack of model energy and to the 
overly large spatial scales used in the Lagrangian analyses . 
The results of the uniform and non-uniform sampling 
strategies used here show good agreement between 
dif fusivities, length, and time scales; but this is 
misleading, since the latitude bands used to group the mean 
drifter positions for statistical analysis cover such a 
large, non-localized meridional scale and span different 
current regimes . 

The low variability, longer spatial scales, low 
dif fusivities, and longer time scales relative to observed 
values indicate that the model fields are missing high 
frequency processes. This is certainly due in part to the 
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form of the wind forcing (3 -day averages of daily ECMWF 
fields interpolated onto the model grid from a 2.5° x 2.5° 
grid) and the use of the 9-day running average in addition to 
any real lack of mesoscale activity in the model. 

B . RECOMMENDATIONS 

The most prominent model weaknesses are weak 
variability, as much as half of that in the observed data and 
excessive mean flows. This is a result of a combination of 
no surface mixed layer, lack of real heat and freshwater 
fluxes, too coarse of a grid resolution, and possibly a 
result of inadequate representation of wind forcing. The use 
of weighted model means, proportional to the number of 
drifter observations per latitude band, would remove some 
ambiguity in the interpretation of differences between model- 
and drifter-based statistics. Incorporation of a surface 
mixed layer would allow the model to achieve an exchange of 
momentum, heat, and freshwater between the ocean and 
atmosphere (represented here by the surface forcing) , and is 
likely to be necessary to correctly simulate the observed 
seasonal cycles in velocity and SST (Stammer et al . , 1996). 
(The most current POP version has implemented a mixed layer 
formulation known as the "K-prof ile parameterization" (KPP) . ) 

Using real heat flux and evaporation - precipation (E-P) 
values, i. e., ECMWF analysis and re-analysis, would allow 
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better simulation of the variability of thermohaline forcing, 
especially the annual cycle. Maltrud et al . (1997) state 

that more realistic surface currents and heat transport occur 
when the model is driven by heat and freshwater fluxes, 
rather than when it is restored to temperature and salinity 
climatology. 

Finer grid resolution would allow better representation 
of points of separation by currents from coastlines, more 
accurate transports through narrow passages, and better 
realization of the baroclinic / barotropic instabilities 
associated with variability. 

Finer grid resolution has been introduced in the North 
Atlantic POP simulation (M. Maltrud, personal communication) , 
where 1/10° resolution and 37 levels produce very accurate 
means and variabilities. SSH variability from TOPEX/ POSEIDON 
compares very well with the model output from this higher 
resolution grid. 

The eddy heat flux in the covariance of velocity and SST 
fields is important for describing the heat balance and SST 
in the equatorial region. This may be a particularly useful 
diagnostic when evaluating the model's ability to simulate 
climate and climate change. The results obtained here point 
to the cold tongue and the northern instability wave regions 
in the eastern tropical Pacific as important locations upon 
which to focus heat budget calculations . A detailed heat 
budget of this region is needed to understand and balance the 
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equatorial cooling, induced by the divergent Ekman transport, 
against the equatorward heat flux by mesoscale eddies . 

As is the case for the Eulerian comparisons described in 
this study, the Lagrangian statistics computed from model- 
derived trajectories also point to focus areas for further 
comparisons and model development. 

Lagrangian time scales from the model trajectories are 
longer than observed and the dif fusivities are smaller than 
observed. Since different techniques were used to determine 
the time scales in this study (e-folding times) versus 
studies based on the drifter data (zero crossing times), 
systematic differences are expected in these integral 
quantities. However, the expected bias between these two 
methods is for shorter time scales from the e-folding method, 
not longer time scales as this study found. More detailed 
computations of Lagrangian statistics are needed to verify 
these differences. On the model side, shorter trajectories 
with durations on the order of several integral time scales 
(a few months) should be used rather than the two-year long 
trajectories in this study. In this way, more trajectories 
would be available for ensembling and their mean positions 
would be more representative of the regional ensemble to 
which they were assigned. On the data side, combined global 
statistics should be computed for comparison with the model 
results rather than the separate seasonal statistics provide 
by Poulain (personal communication) for this study. 
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A great strength of the Lagrangian methods applied to 
model velocity fields is the ability to control the sampling 
strategy. More should be done to exploit this strength using 
shorter trajectories that will better highlight the 
differences between uniform and non-uniform deployments. 

This sampling flexibility is also ideal for the investigation 
of the effects of shear in the mean flow on the computation 
of Lagrangian quantities like integral time scales and 
difusivities (Bauer et al.,1997) as well as Eulerian 
properties like velocity and divergence (Poulain, 1993) and 
mixed sampling problems like the array bias of mean current 
estimates due to non-uniform concentrations in the presence 
of diffusion (Davis et al . , 1996). 
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APPENDIX. 



RUNGE-KUTTA METHOD OF TRAJECTORY GENERATION 



In this study we consider an initial value problem which 
consists of a first-order differential equation and a 
condition which the solution must satisfy (Nakamura 1991; 
Kreyszig 1986) : 

y'=f(y,t), y(0) = y 0 . 

The estimate y n+1 at £ fl+1 = t n + h with a known value of y n , may 
be obtained by integrating over the interval [t n , t ]: 

y n+ \ = yn+ U(y>t)dt. 

‘n 



Following Press et al. (1992): 

A practical method for solving initial value 
problems is the Runge-Kutta (RK) method, which 
propagates a solution over an interval by combining 
the information from several Euler-style steps 
(each involving one evaluation of the right hand 
f's) , and then using the information to match a 
Taylor's series expansion up to some higher order. 

The formula for the Euler method is 

y n + l =yn + h f( x n’y n't 

which advances a solution from x n to x n+1 = x n+h . 

The practical value of the Euler method is limited, 
however, since it advances a solution through an 
interval h, but uses derivative information only at 
the beginning of that interval (Figure 48) . The 
step's error, therefore, is only one power smaller 
than the correction, and the method is not as 
accurate as other methods, or as stable as is 
usually desired for practical use. 

However, the Euler method is useful to take a 
trial step to the midpoint of an interval. Then 
the value of both x and y at that midpoint is used 
to compute the step across the interval (Figure 
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49) . This symmetrization cancels out the first- 
order error term; the method is called the midpoint 
or second-order Runge-Kutta method. 

In RK methods, the order of accuracy is 
increased by using intermediate points in each step 
interval h. Higher accuracy also implies that 
errors decrease more quickly than in lower-order 
accuracy methods when the step interval is reduced. 
The fourth-order RK method requires four 
evaluations of the right-hand side per step h 
(Figure 50) , and is accurate to the fourth-order 
term of the Taylor expansion, so the local error is 
proportional to h 5, 



For this specific problem, the implementation of the RK 
method using the simulation velocity field involves 

X new = X old + uAt ’ y n ew=yold + V A At = 6 ^ • 



Using 



k l =h u (x, y, t), l l =hv(x,y,t ) 



k l h k l 

k^ = h u (x n +~ 2 ~, y n +~2 , t n + ~ 2 ^ ^2 = k v ( x n + ~ 2 ~ ’ y>i + ~2 ’ 



k 1 u 

2 2 h 



k l 

2 2 h 



k 3 = hu(x n +—,y n + ~, t n +~), l 3 = hv(x n + —,y n + — , t n +~) 



k 4 = hu(x n + k 3 ,y n + l 3 , t n + h), l A = hv (x n + k 3 , y n + l 3 , t n + h) 



the two-dimensional trajectories of the simulated drifters 
are computed using the fourth-order RK scheme, where x^ +1 

represents the new location of a particle that is advected 
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from a previous position x , by the velocity u or v, in a time 
interval At, and k ± and 1 £ represents the RK coefficients. 
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■Vi 



Figure 48. 'Euler's method. In this simplest (and least 
accurate) method for integrating an ODE, the derivative at 
the starting point of each interval is extrapolated to the 
next function value. The method has first-order accuracy. 
Figure 16.1.1 from Press et al., 1996. 
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v(x) 




Figure 49. Midpoint method. Second-order accuracy is 
obtained by using the initial derivative at each step to find 
a point halfwauy across the interval, then using the midpoint 
derivative across the full width of the interval. In the 
figure, filled dots represent final function values, while 
open dots represent function values that are discarded once 
their derivatives have been calculated and used. Figure 
16.1.2 from Press et al . , 1996. 
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Figure 50. Fourth-order Runge-Kutta method. In each step 
the derivative is evaluated four times: once at the initial 
point, twice at trial midpoints, and once at a trial 
endpoint. From these derivatives the final function value 
(shown as a filled dot) is calculated. Figure 16.1.3 from 
Press et al., 1996. 
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